Lattice-Boltzmann Simulations of Fluid Flows in MEMS 
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The lattice Boltzmann model is a simplified kinetic method based on the particle distribution 
function. We use this method to simulate problems in MEMS, in which the velocity slip near the 
wall plays an important role. It is demonstrated that the lattice Boltzmann method can capture 
the fundamental behavior in micro-channel flow, including velocity slip and nonlinear pressure drop 
along the channel. The Knudsen number dependence of the position of the vortex center and the 
pressure contour in micro-cavity flows is also demonstrated. 
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The development of technologies in Micro-electro- 
mechanical systems (MEMS) has motivated the study of 
fluid flows in devices with micro-scale geometries, such as 
micro-channel and micro-cavity flows W . fn these flows, 
the molecular mean free path of fluid molecules could 
be the same order as the typical geometric length of the 
device; then the continuum hypothesis which is the fun- 
damental for the Navier-Stokes equation breaks down. 
An important feature in these flows is the emergence 
of a slip velocity at the flow boundary, which strongly 
affects the mass and heat transfer in the system. In 
micro-channel experiments, it has been observed that the 
measured mass flow rate is higher than that based on a 
non-slip boundary condition @J. The Knudsen number, 
K n = l/L, can be used to identify the influence of the 
effects of the mean free path on these flows, where I is the 
mean free path of molecules and L is the typical length of 
the flow domain. It has been pointed out that for a sys- 
tem with K n < 0.001, the fluid flow can be treated as con- 
tinuum. For K n > 10 the system can be considered as a 
free-molecular flow. The fluid flow for 0.001 < K„ < 10, 
which often appears in the MEMS 0], can not be con- 
sidered as a continuum nor a free-molecular flow. Tradi- 
tional kinetic methods, such as molecular dynamics sim- 
ulations H and the continuum Boltzmann equation ap- 
proach, could be used to describe these flows. But these 
methods are more complicated than schemes usually used 
for continuum hydrodynamic equations. The solution 
of the Navier-Stokes equation including the velocity-slip 
boundary condition with a variable parameter has also 
been used to simulate micro-channel flows jl| . 

In the past ten years, the lattice Boltzmann method 
(LBM) j5j has emerged as an alternative numerical tech- 
nique for simulating fluid flows. This method solves a 
simplified Boltzmann equation on a discretized lattice. 
The solution of the lattice Boltzmann equation converges 
to the Navier-Stokes solution in the continuum limit 
(small Knudsen number). In addition, since the lattice 
Boltzmann method is intrinsically kinetic, it can be also 
used to simulate fluid flows with high Knudsen numbers, 
including fluid flows in very small MEMS. 

To demonstrate the utility of the LBM, we use the 
LBM model with three speeds and nine velocities on a 
two-dimensional square lattice. The velocities, Ci, in- 
clude eight moving velocities along the links of the square 



lattice and a zero velocity for the rest particle. They are: 
(±1,0),(0,±1),(±1,±1),(0,0). Let /i(x,i) be the dis- 
tribution functions at x, £ with velocity Ci. The lattice 
Boltzmann equation with the BGK collision approxima- 
tion |(||7j can be written as 

/i (x,f) = -r- 1 (/ i -/n, (1) 



/i(x + Cifa,t + 8t) 



where f^ q (i — 0, 1, • • -,8) is the equilibrium distribution 
function and r is the relaxation time. We have assumed 
that the spatial separation of the lattice is Sx and the 
time step is St. A suitable equilibrium distribution is [Q: 
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Here c s = 1/V3, t Q = 4/9, ti = t 2 = £3 = £ 4 = 1/9 
and £5 = te = £7 = t% = 1/36. The Greek subscripts a 
and (3 denote the spatial directions in Cartesian coordi- 
nates. The density p and the fluid velocity v are defined 
by p — Ylifii P w = J2i c ifi- I n previous lattice-BGK 
models, r was chosen to be a constant. This is applica- 
ble only for nearly-incompressible fluids. In micro-flows, 
the local density variation is still relatively small, but the 
total density change, for instance the density difference 
between the inlet and exit of a very long channel, could 
be quite large. To include the dependence of viscosity on 
density we replace r in Eq.(|l|) by r': r' = | + jj(r — \). 
Using the Chapman-Enskog multi-scale expansion tech- 
nique, we obtain the following Navier-Stokes equations 
in the limit of long wavelength and low Mach number: 



d t p + d a (pu a ) = 0, 
d t {pu a ) + dp(pu a up) = -d a P 



(3) 
(4) 



P = c Ip, ^ap = v{d a {pup) + dp(pu a )), 



where v = c 2 s (2t — l)/(2p) is the kinematic viscosity. In 
classical kinetic theory, the viscosity v for a hard sphere 
gas is linearly proportional to the mean free path. Sim- 
ilarly, we define the mean free path I in the LBM as: 
a(r — 0.5) / p, where a is constant. 

Our first numerical example is a micro-channel flow [Q . 
The flow is contained between two parallel plates sepa- 
rated by a distance H and driven by the pressure differ- 
ence between the inlet pressure, Pi, and exit pressure, P e . 
The channel length in the longitudinal direction is L. We 
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take L — 1000, H = 10 (lattice units) in our simulations 
satisfying L/H » 1. The bounce-back boundary condi- 
tion is used for the particle distribution functions at the 
top and bottom plates, i.e., when a particle distribution 
hits a wall node, the particle distribution scatters back 
to the fluid node opposite to its incoming direction. A 
pressure boundary condition is used at the input and the 
exit. 
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FIG. 1. The slip velocity and the normalized mass flow 
rate at the exit of a micro-channel flow as functions of K n for 
Pi/Pe = 2. The ' + ' and 'x' are LBM numerical results. The 
dashed and dotted lines are Eq.(^) and Eq.(^) respectively. 

The slip velocity V s at the exit of the micro-channel 
flow is defined as: u(y) — u (Y — Y 2 + V s ), where u(y) is 
the velocity along the x (or flow) direction at the exit and 
Y = y/H. uo and V s can be obtained by fitting numerical 
results using the least squares method. This definition of 
the slip velocity is consistent with others (|,||]. I n Fig- 1, 
we plot the slip velocity V s and the normalized mass flow 
rate Mf = M/M , as functions of Knudsen number when 
the pressure ratio rj = Pi/P e = 2. The normalization fac- 
tor, Mo — ^arf fa — 1); is the mass flow rate when the 
velocity slip is zero. To calculate the Knudsen number 
we have chosen a = 0.388 in order to match the simu- 
lated mass flow rate with experiments (See theory curve 
in Fig. 3). Using a least squares fit to the data in Fig. 1, 
we obtain: 



V s = %.!Ki 



(5) 



If we assume that the Navier-Stokes equations are valid 
for the micro-flows except that the slip boundary condi- 
tion V s in Eq.^) replaces the traditional non-slip condi- 
tion on the walls |j, Eqs. (§), (§) and (!) wil1 S ive the 
mass flow rate: 



Mf = 1 + 12V S (K„ 
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(6) 



For j] — 2, the above formula becomes Mf = 1 + 24.1K 2 , 
which agrees well with the numerical results in Fig 1. 

In laminar Poiseuille flows, one usually assumes that 
the density variation along the channel is very small, and 
the pressure drop along the channel is nearly linear. In 



micro-channel flow, however, the ratio between the length 
and the width is much larger and the pressure drop is 
not linear. If there is no velocity slip at the walls, it has 
been shown from the Navier-Stokes equation that 
the pressure along the channel has the following depen- 
dence on the dimensionless coordinate, X = x/L: 



P 



P e 2 [l + (r, 2 -l)(l-X)], 



(7) 



If the velocity at the boundaries is allowed to slip, the 
pressure drop along the channel will depend on the Knud- 
sen number. In Fig. 2 we present the LBM simulation re- 
sults for the normalized pressure deviation from a linear 
pressure drop, (P — P{)/P e , as functions of X for several 
Knudsen numbers, where Pi = P e + (Pi — P e )(l — X). 
It is seen that when K n < 0.2, (P — P{)/P e is a positive 
nonlinear function of X. This agrees with the results 
in [Q using an engineering model. For K n > 0.2, the 
LBM simulation shows that (P — P{)/P e becomes neg- 
ative, which is directly linked to the fact that the slip 
velocity depends on the square of K n in the LBM. For 
large K n , the pressure can be derived from Eq.(0): 



P = P e [V {1 - X) ] 



(8) 



The negative deviation from a linear pressure drop has 
not been experimentally observed before and it would be 
interesting to testify this experimentally. 
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FIG. 2. The deviations from linear pressure drop for 
■q — Pi/Pe — 2. The top and bottom lines are the analyt- 
ical results from Eq.(g) for K n = and Eq.Q for K „ » 1 
respectively. The other curves are LBM numerical results for 
the Knudsen numbers indicated. 

In Fig. 3 the mass flow rates as functions of the pres- 
sure ratio r\ when K n = 0.165 are shown for our theory, 
the experimental work Q , the engineering model [Q] and 
the LBM simulation. Our theory and the LBM simula- 
tion agree well with the experimental measurements. It is 
noted that for large pressure ratios (r) > 1.8), the LBM 
agrees reasonably well with Beskok et al. Q. But for 
smaller pressure ratios, the difference increases because 
of different dependence of the slip velocity on K n . 
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FIG. 3. The normalized mass flow rate as a function of the 
pressure ratio for K n = 0.165. The theory is Eq.M. 
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FIG. 4. Streamlines for two Knudsen numbers. 
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FIG. 5. The y-coordinate of the vortex center (square sym- 
bols) and the mass flow(solid circles) as a function of K n . 

Our second LBM numerical simulation is the two- 
dimensional micro-cavity flow. The cavity size is L x = 
L y = 40 (lattice units). The upper wall moves with a 
constant velocity, i>o, from left to right. The other three 
walls are at rest, and bounce-back boundary conditions 
are used. To see the dependence on the Knudsen num- 
ber in our simulations, we fixed the Reynolds number, 
Re = = 2 - 4 x 10 ~ 4 and re q uire the Mach number 
to be small, M a = < 1(T 3 . In Fig. 4, wc show the 
streamlines for two different Knudsen numbers. In Fig. 5 
we show the vertical positions of the vortex center and 
the mass flux between the bottom and the vortex cen- 



ter as functions of K n . It can be seen that the center 
moves upward and the mass flow decreases with increas- 
ing Knudsen number. This occurs because the slip ve- 
locity on the upper wall causes momentum transfer to 
be less efficient. It has been shown M that the center of 
the vortex moves downward when the Reynolds number 
increases for very small K n , Fig. 6 shows the pressure 
contours for the same parameters as in Fig. 4. Totally 
different pressure structures are observed for these two 
cases. When the Knudsen number is small, the contin- 
uum assumption is valid and the pressure contours are 
almost circles with centers at the left or the right corners. 
On the other hand, due to the slip velocity on the walls, 
the pressure contours become nearly straight lines at the 
higher Knudsen number. 
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FIG. 6. The contours of pressure for the same two Knudsen 
numbers shown in Fig. 4. 
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